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We study the resistance fluctuation distribution of a thin film near electrical breakdown. The 
film is modeled as a stationary resistor network under biased percolation. Depending on the value 
of the external current, on the system sizes and on the level of internal disorder, the fluctuation 
distribution can exhibit a non-Gaussian behavior. We analyze this non-Gaussianity in terms of the 
generalized Gumbel distribution recently introduced in the context of highly correlated systems near 
criticality. We find that when the average fraction of defects approaches the random percolation 
threshold, the resistance fluctuation distribution is well described by the universal behavior of the 
Bramwell-Holdsworth-Pinton distribution. 



I. INTRODUCTION AND MODEL 

Electrical breakdown is an irreversible change in the re- 
sistance of a two terminal device and it usually occurs in 
the presence of high applied electric stress, like current or 
voltage. Breakdown of condunctor-insulator composites 
[1], soft dielectric breakdown in ultrathin oxide films [2] 
and degradation of metallic films due to electromigration 
[2,3] are typical experimental examples. Nonlinear cur- 
rent voltage characteristics, giant enhancement of excess 
noise over the quadratic dependence on the applied field 
are the mostly studied breakdown precursors. [1,4]. Yet 
another important and less understood breakdown pre- 
cursor is represented by the emergence of non-Gaussian 
fluctuations [5] , which origin has attracted an increasing 
interest in the recent literature [6-9]. For these reasons, 
here we study the distribution of the resistance fluctu- 
ations of a thin film near electrical breakdown. To this 
purpose we make use of the Stationary Network Under 
Biased Percolation (SNUBP) model [3,4]. This model al- 
lows the study of the electrical conduction of disordered 
materials over the full range of applied bias values, from 
the linear regime up to the breakdown. Moreover, the 
SNUBP provides a good modeling of many features as- 
sociated with the electrical instability of composites ma- 
terials [1] and with the electromigration damage of metal 
lines [3], two important classes of breakdown phenomena. 

In the SNUBP model [3,4], the resistance and the resis- 
tance fluctuations of a disordered film are determined by 
the competition of two biased processes taking place in 
a two-dimensional square-lattice resistor network. More 
precisely, the network consists of 2N 2 resistors in two 
possible states: (i) regular, corresponding to resistors 
with resistance r n = r [l + a(T n — T )] and (ii) broken, 
corresponding to resistors with rop = 10 9 r n . Here, N 
determines the linear sizes of the network, a is the tem- 
perature resistance coefficient, To the bath temperature, 



and T„ the local temperature of the n-resistor, resulting 
from Joule heating and thermal exchange with neighbour 
resistors [10]: 
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where N ne i g is the number of first neighbors around the 
nth resistor and i n the current flowing in it. The param- 
eter A represents the thermal resistance of each resistor 
and it determines the importance of Joule heating effects. 
The network is contacted at the left and right hand sides 
to perfectly conducting bars through which a constant 
stress current / is applied. The two biased processes 
consist of stochastic transitions between the two possible 
states and they occur with thermally activated probabil- 
ity, Wd and Wr, which depend on the local temperature 
and are characterized by two activation energies, Ed and 
Er, with Ed > Er for networks of finite sizes [4,11]. The 
network evolution is obtained by a Monte Carlo simula- 
tion which updates the network resistance after breaking 
and recovery processes according to an iterative proce- 
dure as detailed in Ref. [4]. The sequence of the suc- 
cessive configurations provides the final R(t) signal with 
the time scale calibrated on the iteration steps. Then, 
depending on the Ed and Ep values and on the stress 
conditions, the network reaches a steady state or under- 
goes an irreversible electrical failure. This last possibility 
is associated with the existence of at least one continu- 
ous path of broken resistors between the upper and lower 
sides of the network, i.e. with the achievement of the 
percolation threshold, pc, for the fraction of broken re- 
sistors. It has been found [4] that for a given value of To, 
a threshold current value, Ib , exists above which the elec- 
trical breakdown occurs. On the other hand, for I < Ib 
the network is stable, i.e. the fraction of broken resistors, 
p, and the network resistance, i?, fluctuate around their 
average values < p > and < R >, respectively. In the 
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following section we will analyze the results of simula- 
tions carried out by considering square networks of size 
75 x 75 and by taking the following values for the other 
parameters (chosen as physically plausible): ro = 1 (Q), 
a = 1(T 3 (K- 1 ), A = 5 x' 10 5 (K/W), E D = 0.170 (eV), 
E R in the range 0.026^-0.155 (eV) and T = 300 (K). The 
values of the external current range from 0.001 < I < 3.0 
(A). 

II. RESULTS 



Figure 1 reports a set of resistance evolutions for in- 
creasing current values. In this case the activation energy 
for the recovery process is Er = 0.043 (eV). This value 
of Er implies a breakdown current Ib = 2.1 (A). The 
figure displays two important features of the electrical 
response of conducting films to high fields: i) the aver- 
age resistance exhibits an Omic behavior at low currents 
and then increases significantly at high currents (non- 
Ohmic regime); ii) the amplitude of resistance fluctua- 
tions increases strongly for currents near to the break- 
down value. By analyzing the first and the second mo- 
ment of the steady resistance we have found that the 
non ohmic regime can be divided in two regions [4] . Ac- 
tually, a moderate bias region, where both the average 
resistance < R > and the relative variance of resistance 
fluctuations S =< (Ai?) 2 > / < R > 2 increase quadrat- 
ically with the current, is followed by a pre-breakdown 
region where superquadratic behaviors emerge [4] . 
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FIG. 1. Resistance evolutions for increasing bias values. 
Starting from the bottom, the six curves correspond to steady 
states at I = 0.01, 0.5, 1.0, 1.5, 1.8, 2.1 (A). The highest curve 
corresponds to breakdown and it is obtained for / = 2.2 (A). 

Remarkably, another feature of the resistance fluctua- 
tion distribution is found to emerge: the non-Gaussian 
behavior which accounts for the importance of moments 
higher than the second. Here, we have thus investigated 
the Gaussianity of the R(t) steady-state signals shown in 
Fig. 1. We have considered time series containing about 
3 x 10 5 resistance values. The results of the analysis are 
reported in Fig. 2 which shows the probability density 
function (PDF), </>, of the distribution of resistance fluc- 
tuations for two values of the bias current: / = 0.01 (A) 
(triangles) and I = Ib = 2.1 (A) (full circles). Precisely, 



in this figure we have reported on a lin-log plot the prod- 
uct cr<I> as a function of (< R > —R)/<r (where er is the 
root mean square deviation). In fact, by making the dis- 
tribution independent of its first and second moments, 
this normalized representation is particularly convenient 
to explore the functional form [6]. 
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FIG. 2. Scaled PDF of resistance fluctuations, II = <j(j>, for 
two current values: I = 0.01 (A) (triangles) and I = 2.1 (A) 
(full circles). The other curves are respectively: BHP distri- 
bution (thick solid line), Gaussian (long dashed), best- fit to 
the I = 0.01 (A) data (solid), best-fit to the I — 2.1 (A) data 
(dashed) . 

For comparison, in Fig. 2 we also show the Gaussian 
distribution (long dashed line), which in this normalized 
representation has zero mean and unit variance. We no- 
tice that in 1998 Bramwell, Holdsworth and Pinton [6] 
realized that the PDF of a global quantity of a system at 
a critical point takes a universal behavior, irrespectively 
of the particular quantity and of the nature and sizes 
of the system. Successively, the functional form describ- 
ing this universal behavior was identified [7]. Therefore, 
in Fig. 2 we report also the Bramwell, Holdsworth and 
Pinton (BHP) distribution (thick continuos line). Before 
giving its expression, we have to introduce the following 
definitions [7] . We call m a fluctuating quantity (for ex- 
ample the magnetization of a ferromagnet like in Ref. [7]), 
< m > and a m its mean value and root mean square devi- 
ation respectively, P(m) its PDF, y = (m— < m >)/u m 
the normalized variable, U(y) = a m P(y) the normalized 
PDF and finally x = b(y - s). Then, the BHP distribu- 
tion takes the following expression [7]: 

n(y) = K[e x ~ e *] a (2) 

where a = tt/2, b = 0.936 ± 0.002, s = 0.374 ± 0.001 
and K = 2.15 ± 0.01. Equation (2) can be considered a 
generalization of the Gumbcl distribution, which is often 
associated with the occurrence of rare events. Thus, look- 
ing at Fig. 2 we can outline the following main points: 
i) the resistance fluctuation distribution exhibits a non- 
Gaussian behavior; ii) the non-Gaussianity increases for 
currents near to the breakdown value; iii) even at low bi- 
ases (i.e within the Ohmic regime) the distribution keeps 
a weak non-Gaussianity; iv) the distribution at Ib de- 
parts significantly from the universal behavior. Point iii) 
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is understood as a finite size effect, as confirmed by the 
analysis of data for networks of different sizes, not re- 
ported here. Thus, the points iii) and iv), by indicating 
a size dependent and non universal behavior, both wit- 
ness a first order transition, in agreement with the results 
in Refs. [4,12] and with electrical breakdown experiments 
in composites [1]. On the other hand, by recalling that 
real systems have finite sizes, the feature ii) can provide 
an important precursor of failure for systems of finite and 
given sizes. In order to identify what kind of distribution 
describes the data in Fig. 2, we have considered several 
distributions among those available in the literature. Ac- 
cordingly, we have found that the best description to the 
present data is given by the expression in Eq. (2) once 
the parameters a, b, s and K are considered as fitting 
parameters. The dashed and the solid curves in Fig. 2 
represent the best- fit with Eq. (2) to the data for J = 2.1 
(A) (full-circles) and for / = 0.01 (A) (triangles), respec- 
tively. 
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FIG. 3. Scaled PDF of resistance fluctuations, II = o~4>, for 
different values of the recovery energy. The full circles corre- 
spond to Er = 0.043 (eV), the stars to Er = 0.130 (eV) and 
the squares to Er = 0.155 (eV). In all cases I = Ib- The 
other curves are respectively: BHP distribution (thick solid 
line), Gaussian (long dashed), best-fit to the Er — 0.043 (eV) 
data (dashed), best-fit to the Er = 0.155 (eV) data (solid). 

In the following, we further consider the effect of the 
activation energy of the defect recovery process on the 
PDF. In fact, the recovery energy, by controlling the level 
of disorder inside the system, plays an important role on 
the system behavior [4]. To this purpose, Fig. 3 reports 
the distributions of the resistance fluctuations obtained 
for different values of the recovery energy: E R = 0.043 
(eV) (full circles), E R = 0.130 (eV) (stars), E R = 0.155 
(eV) squares. In all the cases, the current corresponds to 
the breakdown, and the values are: I = 2.1 (A), I = 0.38 
(A) and I = 0.095 (A), respectively. Thus, full circles are 
the same data of Fig. 2. We notice that E R = 0.155 (eV) 
is very close to the Ed value and practically represents 
one the largest value of E R yet providing a steady state 
of the network [4]. Therefore, as discussed in Ref. [4], 
for this value of E R the network is very near to its crit- 
ical point and it is characterized by an average fraction 
of defects (< p >= 0.35) close to the value of the ran- 



dom percolation threshold (0.5 for bond percolation on a 
square and infinite lattice). Accordingly, from Fig. 3 we 
can see that the resistance fluctuation distribution cor- 
responding to this value of E R is well described by the 
universal BHP distribution. The solid and the dashed 
curves in this figure represent the best-fit with Eq. (2) to 
the data obtained for E R = 0.130 (eV) and E R = 0.043 
(eV), respectively. Also in these cases the fit with Eq. (2) 
is found to be satisfactory. 

In conclusions, we have studied the distribution of re- 
sistance fluctuations of a thin film under different bias 
conditions and for several values of the recovery energy. 
The distribution exhibits a non-Gaussian behavior which 
is well described by a generalized form of the Gumbel dis- 
tribution. Furthermore, a value of the recovery energy 
close to that of defect generation, a condition which cor- 
responds to an average defect fraction approaching the 
random percolation value, is identified as the condition 
under which the distribution achieves the universal be- 
havior described by the BHP distribution [6]. 
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